clear
clc

% 参数命名规则：(temp)_数据意义_数据类型
% 数据类型：Tspace 间隔时间
%          Series 其他序列
%          Martrix 矩阵
%          Tseries 时间序列
%          Freq 频率
%          Num 数量
%          Cnt/Sum 计数

sign_Tspace = 1; % 符号间隔
carrier_Freq = 5 / sign_Tspace; % 载波频率，一个符号里面有40个载波周期
sample_Tspace = 1 / (5 * carrier_Freq); % 采样间隔，表明一个周期里面采样5次
sample_Freq = 5 * carrier_Freq; % 采样频率
bit_Num = 2048; % 传输8192个符号
bit_Series = randi(2, bit_Num, 1)' - 1; % 传输的比特序列
sample_Tseries = 0 : sample_Tspace : sign_Tspace * bit_Num; % 采样的时间点序列

modulate_Level = 1; % 调制水平为1，表明0信号的调制幅度为0
SNR_Series = 10.^((-8: 20) ./ 10); % 信噪比幅度表示序列 

% 载波
modulate_AMPspace = 1; % ASK的调制间隔幅度
ASK_Output = zeros(1, length(sample_Tseries)); % ASK的输出结果
X = bit_Series * modulate_AMPspace; % ASK映射后的每点幅度
gray_MapTable = [0;1]; % 格雷码映射表

% 非相干解调滤波器
FIR_Filter = [
  -0.000449919922023
  -0.010290061250890
  -0.020740302584492
  -0.021666933945834
   0.004057013454896
   0.063826480507302
   0.144878533787183
   0.216082197967452
   0.244509350528873
   0.216082197967452
   0.144878533787183
   0.063826480507302
   0.004057013454896
  -0.021666933945834
  -0.020740302584492
  -0.010290061250890
  -0.000449919922023
];

for cnt = 1: bit_Num
    index_Series = ((cnt - 1) * sample_Freq + 1) : (cnt * sample_Freq); % 对一个符号周期内的点输出
    temp_g(index_Series) = 1; % 加权幅度
    % 针对每个符号对应的序列，生成ASK输出
    ASK_Output(index_Series) = X(cnt) .* temp_g(index_Series) .* cos(2*pi * carrier_Freq .* sample_Tseries(index_Series));
end
plot(sample_Tseries,ASK_Output);%绘制调制后图像
xlim([0,10]);
ylim([-2,2]);

ASK_Half = ASK_Output;
ASK_Half(ASK_Output < 0) = 0;
plot(sample_Tseries, filter(FIR_Filter, 1, ASK_Half));
xlim([0,10]);
ylim([-2,2]);

% 添加噪声
SER = zeros(1, length(SNR_Series));
BER = zeros(1, length(SNR_Series));
MC_Times = 10; % 执行100次蒙特卡洛模拟
Noise_delta_Series = 0.5 * modulate_AMPspace ^ 2 ./ SNR_Series; 
X_out = zeros(1, length(bit_Num)); % 执行了AWGN判决以后的输出结果

for cnt = 1:length(SNR_Series)
    fprintf("测试信噪比%f\n", 10 * log10(SNR_Series(cnt)));
    for cnt_mc = 1:MC_Times
        Noise_Series = sqrt(Noise_delta_Series(cnt)/2) .* (randn(1, length(ASK_Output)) + 1i * randn(1, length(ASK_Output))); % 噪声时域函数
        
        ASK_Noise_Out = ASK_Output + real(Noise_Series); % 就是混合了噪声以后的输出序列，对应于X_out是从中判决生成的结果数据
        
        figure(cnt * 4)
        plot(sample_Tseries,ASK_Noise_Out);%绘制调制后图像
        xlim([0,10]);
        ylim([-2,2]);
        
        % 在这里执行半波整流并且低通滤波
        ASK_Noise_Out(ASK_Noise_Out < 0)=0;
        
        ASK_Filter_Out = filter(FIR_Filter, 1, ASK_Noise_Out);
        
        figure(cnt * 4 + 3);
        plot(sample_Tseries,ASK_Filter_Out);%绘制调制后图像
        xlim([0,10]);
        ylim([-2,2]);
        
        for cnt_Y = 1:bit_Num
            index_Series = ((cnt_Y - 1) * sample_Freq + 1) : (cnt_Y * sample_Freq);
            temp_sum = length(find(ASK_Filter_Out(index_Series) > 0.15));
            if (temp_sum >= sample_Freq * 0.8)
                X_out(cnt_Y) = 1;
            else
                X_out(cnt_Y) = 0;
            end 
        end
        SER(cnt) = SER(cnt) + length(find(X ~= X_out)) / bit_Num;
        BER(cnt) = SER(cnt); % 注意，只有在2ASK的时候才能这么干，不然还要加上解码的东西
    end
    SER(cnt) = SER(cnt) / MC_Times;
    BER(cnt) = BER(cnt) / MC_Times;
end

figure(114514);
plot((-8:20), 10*log10(SER));
xlabel('信噪比/dB')
ylabel('SER/dB')
